function energyKWh = objectiveMainKWh(designVariables)
% This is the main objective funciton file
% It will return the value for whatever we choose to go with. For instance,
% if we go with dollars then it will return the value in dollars.
% INPUTS
    % designVariables: a vector with the value of the design variables in
    % the form [length width numberOfPanels var4 ... varN]
% OUTPUTS
    % output: the cost or energy or whatever we end up using to measure
    
% Last modification by Rodney

% Written collectively by group members Lamar Julien, Rodney Metoyer, and
% Lisa Worthington.
% License information for any code not written by the group is embedded in
% the respective functions. All software written by the authors listed
% above is not intended for release, but may be released under an
% open license only.

%% Simulation Set-up

% Parse the design vector
fieldLength = designVariables(1);
fieldWidth  = designVariables(2);
panelLength = designVariables(3);
panelWidth  = designVariables(4);
panelSpace  = designVariables(5);

% Set the simulation constants
rowSpace   = 1.0;   % The north/west spacing between rows
dayID      = 1;     % The (Julien?) day number
latitude   = 35.78; % The latitude for the simulation (Raleigh)
gravity    = 9.81;   %meters/sec^2
panelThick = 0.0005; %meters
densityAir = 1.2133; %Standard atmosphere Raleigh
densPanel  = 24454.5; % Denstity of monocrystalline silicone plus 5% for framing

% Initialize sim time
% Simulation begins at solar noon because we have symmetry (sunrise to noon
% and noon to sunset)
% NOTE: time is in hours and all time is expressed in "solar time" with
% true zenith defined as 12.
simTime = 12.0;
aSecond_hrs = 1/3600;
% aMinute = 1/60;
simDt_hrs = aSecond_hrs;
simDt_sec = simDt_hrs*3600;

% fastDt = aSecond;

% Initialize panel elevation angle (phi)
elevAngle = 0; % In degrees. 0 = flat on the ground
% The torque that the motor is putting out is initially zero, then the
% controller moves it to where it needs to be each time step. We are
% assuming that this particular motor has an electronic lock so that it
% holds at the end of each cycle.
nextElevAngle = 0;
angularVelocity = 0; % Will be estimated
motorTorque = 0;
motorPower = 0;
motorEnergy = 0;
phi = [elevAngle, angularVelocity];
generatedEnergy = 0;

% Initialize wind
wind = struct('direction',1,'speed',1);

% Init misc
flag = true;
sunSetTime = sunset(dayID, latitude);

%% Run Simulation
while(flag)
   % Power Generated by the panel for this elevation angle
   [powerGenerated, flag, numPanels] = Power_Gen(fieldWidth,...
       fieldLength, panelLength, panelWidth, rowSpace, panelSpace, dayID,...
       simTime, latitude, elevAngle, simDt_hrs);
   
   % Get the elevation angle for use elsewhere
   elevAngle = opt_beta2(dayID,simTime,latitude); %degs
   % Convert to radians
   elevAngle = elevAngle*pi/180; %rad
   % And the future angle
   nextElevAngle = opt_beta2(dayID,(simTime+simDt_hrs),latitude); %degs
   % Convert to radians
   nextElevAngle = nextElevAngle*pi/180; %rad
   deltaAngle = nextElevAngle - elevAngle;
   angularVelocity = deltaAngle/simDt_sec; %rad/sec

   % Put it in an array to use later
   phi = [elevAngle, angularVelocity];
   
   % Get the wind
   wind = genWind(wind);
   windDirCode = wind.direction;
   windDir = windDirCode*22.5; % In degrees east of north
   windSpeed = wind.speed;
   windProj = -windSpeed*cosd(windDir - 90); % 90 is because the panels always come up east to west (or west to east)
   
   % Get the torque on the panel at thie angle. This is the minimum that an
   % ideal motor would have to overcome to keep the panel moving at a
   % constant angular velocity
   torqGrav = panelLength*panelWidth*panelWidth*panelThick*densPanel*gravity/2*cos(elevAngle);
   if(windProj==0)
       torqWind = 0;
   else
       if elevAngle > pi/4
           windElevAngle = pi/4;
       else
           windElevAngle = elevAngle;
       end
       torqWind = pi*windElevAngle*densityAir*windProj^3/abs(windProj)*panelWidth*panelLength/4;
   end
   
   % Multiply the torque by the angular velocity estimate to get power
   motorPower = (torqGrav + torqWind)*angularVelocity;
   % Multiplied by the number of panels
   motorPower = motorPower*numPanels;
   
   % Integrate power to get energy
   motorEnergy = motorEnergy + motorPower*simDt_sec;
   generatedEnergy = generatedEnergy + powerGenerated*simDt_sec;
   
   % Advance time
   simTime = simTime + simDt_hrs;
   
   % uh-oh bail out
   if simTime > 100
       disp('We have a problem');
       flag = false;
   end
   % Temproary fix to get around the fact that the flag is not working
   if simTime < sunSetTime
       flag = true;
   end
   
end

energyJoules = generatedEnergy - motorEnergy;
energyKWh = energyJoules*2.7777778*10^-7;